#ifndef _CollinsFragmentation_H_
#define _CollinsFragmentation_H_

//
// CollinsFragmentation.h
//
// Created by TB on 11/21/11.
// Copyright 2011 BNL. All rights reserved.
//

#include <vector>

#include <TVector3.h>

#include "tppmc/DssFragmentation.h"
#include "tppmc/parameters.h"
#include "tppmc/Particle.h"
#include "tppmc/TF2J.h"

class TH2;

namespace tppmc {

/**
 Applies a spin-dependent modulation to fragmentation.
 Fragmentation is only to pions.
 \remark Which Collins parameterisation is this? Rename appropriately.
 \remark Is inheritance from DssFragmentation appropriate? Having a
 member DssFragmentation for unpolarised behaviour may be more natural...
 \todo Need code to handle unsupported parton flavours, namely
 random phi angle.
 */
class CollinsFragmentation : public DssFragmentation {
   //   class CollinsFragmentation : public FragmentationFunction {
   
public:
   
   /**
    Constructor.
    The Collins function applies a transverse-spin-dependent modulation
    to fragmentation.
    \arg \c ff: The fragmentation function which provides the unpolarised
    fragmentation behaviour.
    \arg \c spin: The proton beam spin direction in the lab frame. The
    vector is normalised if the input is not.
    \c ff is not owned by the CollinsFragmentation and so must point to an
    object with lifetime greater than the CollinsFragmentation.
    \todo ff, spin arguments aren't used - either implement functionality
    based on them or scrap them.
    */
   CollinsFragmentation(const FragmentationFunction* ff = NULL,
                        TVector3 spin = TVector3(0., 1., 0.),
                        int order = LO);
   
   
   /**
    Destructor.
    */
   virtual ~CollinsFragmentation();
   
   
   /**
    Fragment the input parton into a hadron at the requested Q2 scale.
    Returns a vector with a single entry.
    Sets the type, precursor parton and z of the returned hadron.
    Sets the azimuthal angle of the hadron with respect to the spin
    direction via px and py.
    All other fields have their default value.
    */
   virtual std::vector<tppmc::Hadron> Fragment(const pythpart&,
                                               double scale);
   
   virtual double Evaluate(const Parton&, const Hadron&, double scale);
   
   /**
    Perform initialisation operations.
    
    Must be called before anything else.
    */
   virtual bool Initialise();
   
   /**
    Returns true if the fragmentation of the input parton type
    is supported by this function.
    */
   virtual bool Supports(const pythpart&) const;
   
   virtual std::string Name() const;
   
   virtual double MinimumScale() const { return 2.; }
   
   /**
    Set either or both the favoured and unfavoured fragmentation functions
    to their positivity bound.
    */
   virtual void SetPositivityBound(bool atBound = true);
   
   TH2* cffPiU();
   
   TH2* cffPi0DbarSd();
   
   virtual void SetHadron(int i = INVALID_HADRON) {
      FragmentationFunction::SetHadron(i);
      mDss.SetHadron(i);
   }
   
   virtual void SetZRange(double lower, double upper) {
      DssFragmentation::SetZRange(lower, upper);
      mDss.SetZRange(lower, upper);
   }
   
   void Qa(std::string filename);
   
protected:
   
   /**
    Implement the hadron list for the collins function as it only
    supports pions, whereas DSS support kaons, protons also.
    */
   virtual const std::vector<int>& Hadrons() const;
   
   
   virtual TF2* Function();
   
   TVector3 RotateToPartonFrame(const pythpart&, const TVector3&) const;
   
   const FragmentationFunction* mUnpolarised;
   TVector3 mSpin;
   PARAMETERS mParameters;
   //      bool mFavouredAtPositivityBound;
   //      bool mUnfavouredAtPositivityBound;
   bool mAtPositivityBound;
   
   DssFragmentation mDss; /// Used to provide unpolarised behaviour
                          /// for unsupported parton types.
   TF2J* mFunction;
   
   ClassDef(tppmc::CollinsFragmentation, 1)
};

inline void CollinsFragmentation::SetPositivityBound(bool atBound) {
   //      mFavouredAtPositivityBound = favoured;
   //      mUnfavouredAtPositivityBound = unfavoured;
   mAtPositivityBound = atBound;
}

} // namespace tppmc


#endif
